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FOREWORD 

This  report  describes  work  performed  in  the  Aerospace  Research 
Laboratories  and  in  the  Air  Force  Flight  Dynamics  Laboratory  under 
the  Defense  Research  Sciences  Program,  Project  7071,  Research  in 
Applied  Mathematics,  Task  01,  Computational  Aspects  of  Fluid  and 
Structural  Mechanics.  The  work  described  was  carried  out  between 
October  1972  and  July  1976.  This  is  an  interim  report.  Further  re- 
ports in  this  series  will  be  documented  under  Project  2304,  Mathe- 
matical and  Information  Sciences,  managed  by  the  Air  Force  Office  of 
Scientific  Research,  Bolling  AFB,  DC. 
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SECTION  I 
INTRODUCTION 

The  structural  analysis  of  mathematical  models  of  complex  aero- 
space and  naval  vehicles  often  involves  the  partial  solution  of  the 
generalized  algebraic  eigenvalue  problem.  A pair  of  n-square  matrices 
A and  B are  given  and  it  is  required  to  calculate  one  or  more  non-zero 
eigenvectors  u and  corresponding  complex  eigenvalues  X for  which  Au  = 
XBu.  In  many  familiar  applications  the  matrices  A and  B are  real  and 
symmetric  and  B is  positive  definite  so  that  all  arithmetic  is  real. 

But  n is  very  large.  Often  A and  B have  relatively  few  non-zero 
entries;  i.e.,  A and  B are  sparse.  Nevertheless,  the  partial  numer- 
ical solution  of  the  large  generalized  symmetric  eigenproblem  is  a 
challenge  to  the  numerical  analyst  and  the  systems  programmer  alike. 

The  terms  "large"  and  "very  large"  are  used  here  relative  to  the  com- 
puter at  hand  and  the  limitations  imposed  by  its  operating  system. 

One  may  distinguish  three  cases.  If  at  most  one  copy  of  A and  B can 
be  stored  in  the  central  memory  of  the  machine,  we  say  that  n is 
large.  It  n is  so  great  that  only  a relatively  few  n-vectors  may  be 
simultaneously  stored,  then  n is  very  large.  If  at  most  one  n-vector 
can  be  accomodated,  then  n is  simply  gigantic . 

This  report  describes  the  simultaneous  iteration  method  for  the 
partial  solution  of  the  generalized  symmetric  eigenproblem  and  gives 
a USA  Standard  FORTRAN  implementation  of  the  algorithm  called  SUBROUTINE 
SIMITZ.  The  use  of  SIMITZ  is  indicated  if  n is  large  or  very  large. 
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hut  it  will  not  handle  gigantic  problems.  However,  it  may 
be  appropriate  to  problems  of  moderate  size  when  only  a few  eigen- 
values or  eigenvectors  are  required. 

Section  II  of  this  report  provides  a brief  mathematical  formu- 
lation of  the  simultaneous  iteration  method  and  describes  the 
Rutishauser-Reinsch  algorithm  ritzit  upon  the  program  SIMITZ  is  based. 
Section  III  is  devoted  to  a discussion  of  alternative  methods  for 
solution  of  the  large  eigenproblem  and  their  comparison  with  simul- 
taneous iteration.  Detailed  instructions  for  use  of  SIMITZ  are  given 
in  Section  IV  and  a complete  listing  of  the  FORTRAN  program  comprises 
Section  V. 
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SECTION  II 
DESCRIPTION 


J 


(.1 


The  present  program  is  an  implementation  of  the  simultaneous 

iteration  algorithm  [2]  for  calculating  the  eigenvalues  largest  in 

magnitude  and  corresponding  eigenvectors  of  a real  matrix  symmetric 

relative  to  a prescribed  inner  product.  Let  ip(n,  w,  z)  denote  an 

inner  product  in  the  space  of  real  column  n-tuples  and  let  the  real 

n-square  matrix  C satisfy  ip(n,  Cw,  z)  = ip(n,  w,  Cz) . Then  C is 

symmetric  relative  to  ip,  and  if  the  n-square  positive  definite 

T 

matrix  B satisfies  ip(n,  w,  z)  = w Bz  then  C is  B-symmetric.  The 
T 

equation  BC  = C B characterizes  the  B-symmetry  of  C.  Given  an  op- 
tional set  of  p initial  approximate  eigenvectors  of  a real  n-square 
B-symmetric  matrix  C corresponding  to  p eigenvalues  of  C largest  in 
magnitude,  the  program  calculates  em  eigenvalues  and  em  corresponding 
eigenvectors,  0 <_ em  < p <.  n,  to  a precision  dependent  on  the  struc- 
ture of  C and  on  a prescribed  tolerance  eps.  The  matrix  B is  pre- 
sented to  the  program  as  an  independently  prepared  real  function  sub- 

T 

program  which  calculates  ip(n,  w,  z)  = w Bz  given  column  n-vectors 
w and  z.  The  matrix  C is  presented  as  an  independently  prepared 
subroutine  subprogram  op(n,  w,  z)  which  when  given  an  n-vector  z com- 
putes its  image  w = Cz.  The  program  is  an  outgrowth  of  a literal 

FORTRAN  translation  [6]  of  the  ALGOL  procedure  ritzit  [9]  to  which 

T 

it  is  substantially  equivalent  when  ip(n,  w,  z)  = w z,  the  standard 
inner  product.  But  depending  on  the  choice  of  B and  C,  the  present 
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program  enables  the  direct  treatment  of  a wide  variety  of  symmetric 
eigenproblems. 

T T 

Let  A = A and  B = B denote  n-square  real  matrices  and  let  a 
be  real.  If  B is  positive  definite  then  the  matrix  C = B 1 (A  - oB) 
is  B-synanetric , and  the  program  computes  eigenvalues  farthest  from 
a of  the  eigenproblem  Au  = XBu  and  corresponding  eigenvectors.  Im- 
plementation of  op(n,  w,  z)  here  consists  in  providing  for  the  appro- 
priate solution  for  w of  the  linear  system  Bw  = (A  - oB)z.  Alter- 
natively, selection  of  op  to  solve  the  system  (A  - aB)w  = Bz  for  w 
enables  the  calculation  by  simultaneous  inverse  iteration  of  the 
eigenvalues  nearest  to  a and  their  eigenvectors.  Implications  for 
large  sparse  systems  for  which  the  Cholesky  factorization  [7]  of  B is 
impractical  are  clear. 

Let  the  eigenvalues  d^,  ...  , d^,  dn+1 , ...  , d^  of  C be  arranged 


in  order  of  descending  absolute  value  and  let  E^  denote  the  direct  sum 


P'  P+1 
ae  and  1 

of  the  distinct  eigenspaces  corresponding  to  dj  , ...  , d^.  Let  Xq 
denote  an  n-by-p  matrix  having  a p-dimensional  column  space  not  orth- 
ogonal relative  to  ip  to  any  eigenvector  in  . Simultaneous  iter- 
ation is  based  on  the  observation  that  if  Id  | > Id  , the  columns 

1 p 1 1 p+1 1 

of  the  matrix  X.  , = C^X,  tend  to  a basis  of  E as  ks  = k + m in- 

k+m  k p 

creases.  But  in  practice  all  of  the  columns  of  X^s  tend  toward  the 
eigenspace  E^  causing  loss  of  information  concerning  the  residual 
eigenvectors.  To  counter  this  tendency,  set 


(1) 
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where  the  p-square  upper  triangular  matrix  B^  Is  constructed  to- 
gether with  Xkg  by  the  Gram-Schmidt  process  to  render  the  columns  of 

X.  orthonormal  relative  to  Ip.  Now  the  i-th  column  vector  of  X, 
ks  r ks 

converges  to  the  i-th  eigenvector  of  C at  a rate  proportional  to 
max  2 < i < p |d^+^/d^|).  Clearly  this  convergence  will 

be  delayed  in  the  presence  of  eigenvalue  clustering.  But  if  |d  | - 
|dp+jJ  is  not  too  small,  the  column  space  of  X^g  will  contain  a good 
approximation  to  the  i-th  eigenvector  even  when  ks  is  small. 

In  order  to  recover  this  approximation,  a modified  Rayleigh-Ritz 

process  is  employed.  Let  denote  an  orthogonal  matrix  which  diag- 

T 

onalizes  the  p-square  symmetric  matrix  bjcsbjcs*  Then  the  i-th  column 


vector  of 


Vi  ' cxkBk«\+i 


(2) 


converges  to  the  i-th  eigenvector  of  C at  a rate  proportional  to 

|dp+^/d^|  while  the  entries  of  the  diagonal  matrix  computed  with  Q^g 

2 2 

and  properly  ordered  offer  close  approximations  to  d, , ...  , d . 

1 P 

The  true  signed  eigenvalues  need  only  be  computed  at  termination  by 

diagonalizing  the  leading  (p  - 1) -square  principal  submatrix  of 
T 

Xj^BCX^,  the  eigenproblem  for  C projected  on  relative  to  ip. 

The  program  determines  a strategy  for  employing  the  devices  (1) 
and  (2)  based  on  the  distribution  of  the  leading  p eigenvalues  of  C 
upon  which  the  convergence  rate  ultimately  depends.  The  selection  of 
values  m in  (1)  is  particularly  important  in  this  regard  in  that  CmX^ 
is  replaced  by  the  m-th  Chebychev  polynomial  on  an  appropriate  inter- 
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val  [-e,  e]  evaluated  by  a special  3-term  recurrence  relation  and 

permitting  accelerated  convergence  when  values  of  m are  continually 

large.  As  a result  the  convergence  quotient  lies  between  |d  /denJ 

and  exp (-arc  cosh  Id  /d  |).  It  is  nearer  to  the  first  value  if 
' em  p 1 

| d ^ / dem | is  large  and  nearer  to  the  second  if  the  latter  quotient 
is  close  to  one. 

As  the  iteration  proceeds  through  a maximum  of  |km|  iteration 
steps  - km  is  a program  parameter  - acceptance  tests  for  the  eigen- 
values and  eigenvectors  are  conducted  following  each  of  the  Rayleigh- 
Ritz  steps  (2).  As  soon  as  the  relative  increase  of  |dji+^|  is  smaller 
than  eps/10,  then  d^+^  is  accepted  and  h,  the  number  of  previously 
accepted  eigenvalues,  is  increased  by  one.  Eigenvectors  are  accepted 
in  groups  of  one  or  more  corresponding  to  clusters  of  accepted  eigen- 
values nearly  equal  in  magnitude.  If  g eigenvectors  have  already 
been  accepted,  let  d , ...  , d.  denote  such  a cluster.  For  all  j, 

g"t"l  X. 

g + 1 j <_  SL,  denote  by  y^  the  projection  relative  to  ip  of  the  image 

Cx.  of  the  j-th  column  x.  of  X,  on  the  linear  closure  of  x,,  ...  , 
j J j Tcs  1 

x . . Set  f . = max  . ||Cx.  — y . | | / | | Cx . | | for  i = g + 1,  ...  , H. 

^ 1 3 J J J 

where  the  indicated  norm  is  the  Euclidean  norm  or  2-norm  relative  to 
ip.  If  |d?  j f J ( |d ^ | - e)  is  smaller  than  eps  then  all  the  x_. , j = 

g + 1,  ...  , Z,  are  accepted  as  eigenvectors  and  g is  increased  to  l. 

The  error  quantities  f^  are  systematically  discounted  in  accordance 
with  the  convergence  properties  of  the  algorithm  to  permit  convergence 
in  the  presence  of  excessive  round-off  error  or  in  case  the  parameter 
eps  is  prescribed  unrealistically  small.  Having  determined  g eigen- 
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vectors,  the  iteration  continues  with  p - g remaining  columns  of 
X^g  until  either  em  eigenvectors  have  been  calculated  or  | km | has 
been  exceeded.  The  program  may  reduce  em  if  it  detects  either  no 
progress  in  convergence  of  eigenvectors  corresponding  to  smaller 
eigenvalues  or  lack  of  stability  in  the  behavior  of  larger  eigen- 
values. 

The  user  may  wish  to  supplement  the  forgoing  outline  of  the 
operation  of  the  program  by  consulting  the  description  of  the  ALGOL 
procedure  ritzit  in  [9]  or  [12]  as  well  as  a review  of  the  mathem- 
atical foundations  of  simultaneous  iteration  in  [5]  and  [8].  For 
this  reason  we  describe  the  principal  differences  between  ritzit  and 
the  present  program.  (a)  The  procedure  inpvod  for  calculating  stand- 
ard inner  products  was  removed  and  the  procedure  ip  was  introduced 
where  appropriate.  (b)  The  procedure  j acobi  for  calculating  the  sol- 
utions of  the  eigenproblem  for  the  p-square  and  (p  - 1) -square  sym- 
metric matrices  was  replaced  by  calls  to  the  EISPACK  [9]  subroutines 
TRED2  and  IMTQL2,  primarily  to  save  space.  (c)  The  procedure  random 
for  calculating  random  column  n-vectors  of  the  matrix  X^g  was  re- 
placed by  in-line  code  which  references  a FORTRAN  function  RANF. 

RANF  returns  uniformly  distributed  random  REAL  values  from  the  inter- 
val (0,  1),  one  per  function  reference,  given  any  one  argument  of 
any  type.  It  is  provided  by  the  user.  (d)  The  procedure  orthog  to 
perform  Gram-Schmidt  orthogonalization  of  the  columns  of  X^g  was  re- 
placed by  internally  linked  in-line  code.  In  attempting  to  control 
potential  underflow  within  orthog  in  a machine  independent  fashion, 
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ritzit  calculates  the  machine  precision  me  but  assumes  in  usage  that 
out-of-range  values  underflow  gracefully  to  zero,  a machine  dependent 
I ' characteristic.  The  present  program  utilizes  a single  REAL  machine 

dependent  constant  MT,  the  ratio  of  the  smallest  FORTRAN  represent- 
able positive  value  to  the  machine  precision,  to  test  for  this  con- 
dition and  upon  its  detection  to  take  appropriate  measures.  (e)  In 

its  ALGOL  implementation  ritzit  requires  approximately  (p  + 3)n  + 

2 

2p  + 5p  storage  locations  in  excess  of  those  required  by  the  program. 

Economies  resulting  in  part  from  (b)  above  [3]  have  reduced  this  re- 

2 

quirement  to  (p  + 2)n  + p + 4p  in  the  present  program.  All  working 

j 2 

storage  is  confined  to  a single  array  of  2n  + p + 3p  locations. 

(f)  The  value  of  km  as  an  input  parameter,  set  to  | km | during  program 
execution,  is  finally  replaced  by  the  value  of  ks  as  an  output  para- 
meter, the  number  of  iteration  steps  used  in  the  calculation  of  em 
eigenvectors.  (g)  The  present  program  retains  unchanged  the  refer- 
ence to  a user  supplied  procedure  inf  as  a window  on  program  execution 
However,  the  one  variable  involving  eps  is  periodically  redefined  to 
enable  effective  control  of  eps  from  inf  or  from  ip  or  op  should  this 
prove  desirable. 

The  testing  procedures  developed  for  the  present  FORTRAN  program 
parallel  its  evolution  from  a research  tool,  which  conformed  closely 
to  its  ALGOL  parent,  to  present  form.  Early  testing  was  concentrated 
on  duplicating  the  tests  furnished  with  ritzit  and  in  eliminating 
errors  in  interpretation  and  translation  of  the  ALGOL  code.  This  was 

: : 

done  for  the  most  part  on  the  CDC  6600  using  FORTRAN  Extended,  Ver- 
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sion  3,  under  the  SCOPE  3.3  operating  system.  Upon  completion  of 
this  first  phase,  the  resulting  program  was  distributed  as  SUBROUTINE 
RITZIT  with  a locally  developed  library  of  FORTRAN  linear  algebra 
routines  [6].  This  same  program  served  as  a basis  of  vitzit  trans- 
lations for  IBM  360/370  processors  [1,3]  whose  preparation  uncovered 
several  bugs  in  the  RITZIT  code  and  suggested  worthwhile  modifica- 
tions. A second  phase  of  testing  involved  the  development  of  a pack- 
age of  auxiliary  FORTRAN  programs  for  use  with  SUBROUTINE  RITZIT  to 
solve  the  eigenvalue  problem  Au  = XBu  through  methods  depending  on 
Cholesky  factorization  where  A and  B may  either  be  full  matrices  or 
sparse  and  banded.  This  phase  was  conducted  on  the  CDC  6600  using 
FORTRAN  Extended,  Version  4,  under  SCOPE  3.4. 

Systematic  testing  of  the  present  program,  SUBROUTINE  SIMITZ, 

has  been  accomplished  in  part  with  the  aid  of  driver  program  TESTB 

which  generates  a symmetric  band  matrix  A and  a lower  triangular 

band  matrix  T of  prescribed  order  n and  bandwidths  whose  relevant 

entries  are  randomly  generated  integer  values  from  a prescribed  in- 

T 

terval.  The  band  matrix  B is  TT  , and  the  program  calculates  the  max- 
imal eigenvalues  of  Au  = XBu.  For  the  sequence  of  values  of  p,  p = 2, 
...  , m in([n/5],  10),  TESTB  exercises  SIMITZ  for  successive  values 
of  em,  em=l,  ...  , p - 1.  For  each  value  of  i,  i = 1,  ...  , em, 
TESTB  computes  the  residuals  Ax^  - d^Bx^  and  their  Euclidean  norms 
relative  to  the  standard  inner  product.  Each  norm  is  normalized  by 
the  difference  |d^|  - e,  and  for  each  value  of  em  the  quantity 

max  , . I I Ax.  - d.Bx. I I /( Id, I - e) , the  value  k of  i for  which 

l<i£emM  l l i 1 1 1 1 1 
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the  maximum  occurs,  and  the  corresponding  geometric  mean  with  unit 
weights  are  listed.  Also  listed  are  the  relevant  non-zero  diagonals 
of  A and  T and  the  final  eigenvalues  computed  for  em  = p - 1. 

Figure  1 shows  an  output  listing  from  the  executable  program 
TESTB  on  the  CDC  6600  under  the  NOS/BE  operating  system  and  FORTRAN 
Extended,  Version  4.  Here  A and  B are  of  order  30  and  each  of  band- 
width 7 having  relevant  entries  between  -99  and  +99.  Listed  are  the 
main  diagonal  and  the  three  adjacent  lower  diagonals  of  T and  A be- 
ginning with  the  entries  in  the  first  column.  Here  eps  = 10  ^ and 
km  = 100.  Note  how  the  relative  nearness  in  magnitude  of  the  first 
three  eigenvalues  inhibits  the  convergence  of  the  second  eigenvector 
when  p = 3 and  em  = 2 resulting  in  acceptance  of  the  first  eigen- 
vector only.  The  fourth  eigenvalue,  however,  is  in  absolute  value 
far  enough  away  from  this  cluster  to  permit  successful  convergence 
when  p = 4 and  em  = 1,  2,  and  3.  This  phenomenon  points  to  a pro- 
cedure for  pursuing  a solution  when  p is  initially  chosen  too  small. 
SIMITZ  may  be  reentered  with  X containing  the  approximate  eigenvectors 
calculated  for  the  smaller  value  of  p as  initial  approximations  for 
use  with  p increased  in  size.  Significant  processor  time  may  often 
be  saved  in  this  fashion. 
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SECTION  III 

SOME  ALTERNATIVE  METHODS 

The  generalized  symmetric  eigenvalue  problem  Au  = XBu  may  in 

theory  be  reduced  to  a standard  symmetric  eigenvalue  problem  Pv  = 

Xv.  When  B is  positive  definite  a lower  triangular  matrix  T may  be 

T 

calculated  by  Cholesky's  method  which  satisfies  B = TT  . Setting 
T -1  -1  T 

v = T u we  have  P * T A (T  ) . This  calculation  of  T may  be 
accomplished  even  when  n is  large.  But  in  practice  the  procedure  is 
complicated  by  the  fact  that  P need  not  be  sparse  even  when  A and  B 
are  so.  But  w = Pz  can  be  calculated  directly  from  A and  T without 
inversion  of  T.  SIMITZ  combined  with  Cholesky's  decomposition  elim- 
inates the  need  for  a special  inner  product  - use  ip  (n,  w,  z)  = 

T 

w z - and  so  improves  efficiency  for  large  n or  for  n of  moderate 
size.  The  so  called  direct  methods  of  treating  Au  = XBu  which  rely 
on  orthogonal  similarity  transformations  to  reduce  P to  tridiagonal 
form  followed  by  application  of  the  QL  or  QR  algorithm  [12]  are  gen- 
erally unsuitable  for  large  problems  unless  the  structure  of  P is 
quite  special.  So  called  iterative  methods  exemplified  by  simul- 
taneous iteration  must  therefore  be  exploited  when  n is  large. 

One  of  the  most  popular  of  such  methods  is  the  inverse  power 
method  or  inverse  iteration  [12].  Given  an  approximate  eigenvalue 
o and  an  approximate  corresponding  eigenvector  the  iteration 

(A  - oB)  = BX.  , , k = 1,  ...  converges  rapidly  to  an  eigenvector 

T 

of  Au  = XBu  when  at  each  stage  X^BX^  = 1.  This  resembles  the  simu- 
taneous  iteration  algorithm  but  when  p = em  = 1 and  when  the  eigen- 
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value  nearest  to  a is  sought  with  its  eigenvector.  The  approximate 
eigenvalue  must  be  calculated  first,  however,  and  inverse  iteration 
may  fail  if  o tends  to  identify  a multiple  eigenvalue.  Eigenvalues 
may  be  obtained  by  any  of  the  various  bisection  methods  [4,12]  or  by 
Rayleigh-Ritz  techniques.  The  inverse  power  method  and  its  variants 
are  the  only  methods  applicable  to  gigantic  problems  known  to  this 
writer. 

One  of  the  most  promising  methods  for  very  large  problems  is 
the  block  Lanczos  algorithm  [11]  for  which  a FORTRAN  program  is  avail 
able  to  solve  the  standard  symmetric  eigenvalue  problem  Pv  = Av. 

Block  Lanczos  has  been  compared  to  ritzit  and  found  to  be  superior  in 
many  situations.  We  see  no  inherent  impediment  to  modifying  block 
Lanczos  to  produce  solutions  of  Au  = XBu  in  the  manner  in  which  we 
have  modified  ritzit. 
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SECTION  IV 

INSTRUCTIONS  FOR  THE  USER 

A FORTRAN  executable  program  or  subprogram  possessing  the  con- 
trol statements  equivalent  to 

DIMENSION  X(MN,P),  D(P),  WK(K) 

INTEGER  P,  EM 

REAL  IP 

EXTERNAL  IP,  INF,  OP 

may  call  SUBROUTINE  SIMITZ  into  execution  via  the  statement 

CALL  SIMITZ  (N,  P,  KM,  EPS,  IP,  OP,  INF,  EM,  X,  MN,  D,  WK) 

where 

N is  an  INTEGER  input  variable,  the  order  of  the  matrix  C. 

P is  an  INTEGER  input  variable,  the  number  of  simultaneous 

iteration  vectors. 

KM  as  an  INTEGER  input  variable  is  in  magnitude  the  maximum 

number  of  iteration  steps  to  be  executed.  If  KM  identifies 
a negative  value  then  P initial  approximate  eigenvectors 
are  assumed  to  be  present  in  the  array  X.  Otherwise  SIMITZ 
supplies  random  initial  eigenvectors. 

KM  as  an  INTEGER  output  variable  identifies  the  number  KS  of 
iteration  steps  finally  used  in  the  calculation  of  EM 
eigenvectors . 

EPS  is  a REAL  input  variable,  the  tolerance  for  accepting 
eigenvectors. 
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IP  is  an  EXTERNAL  input  variable,  the  name  of  a FORTRAN  compat- 
ible REAL  FUNCTION  subprogram  of  the  form  IP(N,  Z,  W)  which 

T 

must  return  the  inner  product  (W,  BZ)  = W BZ  of  the  vectors 
identified  by  the  N-arrays  Z and  W relative  to  the  positive 
definite  matrix  B. 

OP  is  an  EXTERNAL  input  variable,  the  name  of  a FORTRAN  compat- 
ible SUBROUTINE  subprogram  of  the  form  OP(N,  Z,  W)  which  must 
calculate  the  image  W of  the  vector  identified  by  the  N-array 
Z under  the  N-square  matrix  C without  overwriting  Z. 

INF  is  an  EXTERNAL  input  variable,  the  name  of  a FORTRAN  compat- 
ible SUBROUTINE  subprogram  which  may  be  used  for  obtaining  in- 
formation or  to  exert  control  during  execution  of  SIMITZ.  INF 
has  the  form  INF(KS,  G,  H,  F)  where 

KS  is  an  INTEGER  output  variable,  the  number  of  the  next 
iteration  step. 

G is  an  INTEGER  output  variable,  the  number  of  already 
accepted  eigenvectors. 

H is  an  INTEGER  output  variable,  the  number  of  already 
accepted  eigenvalues. 

F is  a REAL  output  variable  P-array,  error  quantities 
measuring  respectively  the  state  of  convergence  of 
the  P simultaneous  iteration  vectors.  In  addition, 
if  convergence  fails  in  SUBROUTINE  IMTQL2  after  G 
eigenvectors  have  been  accepted,  then  F(G+1)  is  re- 
placed by  1000.*FLOAT  (IERR)  where  IERR  is  the  error 
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indicator  output  by  IMTQL2.  Each  element  of  the  array 
F is  initially  set  by  SIMITZ  to  the  value  4.0. 

EM  as  an  INTEGER  input  variable  is  the  number  of  eigenvalues  to  be 
computed,  0 -LT.  EM  .LT.  P .LE.  MN. 

EM  as  an  INTEGER  output  variable  is  the  number  of  eigenvectors  com- 
puted through  KM  iteration  steps. 

X as  a real  N-by-P  input  array  is  a set  of  P optional  initial 

approximate  eigenvectors  X(I,1),  ...  , X(I,P),  1=1,  ...  , N, 
interpreted  by  SIMITZ  if  KM  is  negative. 

X as  a real  N-by-P  output  array  is  a set  of  EM  eigenvectors 

X(I,1),  ...  , X(I,EM),  1=1,  ...  , N,  computed  through  IABS(KM) 

iteration  steps  with  the  remainder  of  X consisting  of  P - EM 

T 

approximate  eigenvectors.  The  N-by-P  matrix  X satisfies  X BX  = 

I,  that  is,  the  eigenvectors  of  C are  B-orthonormal. 

MN  is  an  integer  input  variable  which  identifies  the  leading  dimen- 
sion in  the  calling  program  of  the  array  X. 

D is  a real  output  P-array  of  which  D(l),  ...  , D(EM)  are  the 

eigenvalues  of  C largest  in  magnitude  in  decreasing  order  cor- 
responding to  the  computed  eigenvectors  X(I,1),  ...  , X(I,EM), 
1=1,  ...  , N.  D(EM+1),  ...  , D(P-l)  contain  approximations 
to  progressively  smaller  such  eigenvalues.  D(P)  contains  the 
most  recently  computed  value  of  E,  where  the  interval  (-E,E)  is 

the  interval  over  which  the  Chebyshev  acceleration  was  performed. 

2 

WK  the  initial  location  of  at  least  P + 3P  + 2N  = K consecutive 
storage  locations  which  may  not  be  overwritten  while  SIMITZ  is 
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in  execution. 

SIMITZ  employs  a DATA  statement  to  assign  to  a machine  dependent  REAL 
variable  MT  the  quotient  of  the  smallest  positive  REAL  value  repre- 
sentable by  FORTRAN  and  the  smallest  REAL  value  whose  sum  with  1.0 
exceeds  1.0.  The  performance  of  SIMITZ  is  strongly  dependent  upon 
the  choice  of  input  parameters  and  upon  the  careful  preparation  of 
the  subprograms  IP  and  OP.  The  user  should  develop  experience  with 
SIMITZ  on  problems  of  moderate  size  before  investing  processor  time 
on  very  large  problems  for  which  the  procedure  is  ultimately  intended. 
During  its  execution  SIMITZ  issues  calls  to  the  following  subprograms 
FUNCTION  RANF 

returns  uniformly  distributed  random  numbers  on  the  open 
interval  (0,  1)  given  any  one  argument  of  any  type. 
SUBROUTINE  TRED2 

is  the  EISPACK  program  which  computes  a Householder  tri- 
diagonal form  of  a real  symmetric  matrix. 

SUBROUTINE  IMTQL2 

is  the  EISPACK  program  which  computes  the  eigenvalues  and 
orthonormal  eigenvectors  of  a symmetric  tridiagonal  matrix. 
FUNCTION  IP 

is  described  above. 

SUBROUTINE  OP 

is  described  above. 

SUBROUTINE  INF 

is  described  above. 
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The  Function  RANF  furnished  with  the  CDC  FORTRAN  Common  Library 
Mathematical  Routines  is  an  excellent  random  number  generator.  How- 
ever, no  failures  of  SIMITZ  have  been  noted  even  when  RANF  was  re- 
placed by  a crude  in-line  congruence  procedure.  The  user  should  use 
the  USA  Standard  FORTRAN  versions  of  TRED2  or  TMTQL2  rather  than  the 
DOUBLE. PRECISION  versions  available  for  the  IBM  360/370  compatible 
processors.  The  E1SPACK  subroutine  TQL2  may  replace  IMTQL2  if  de- 
sired. 
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SECTION  V 

« ALGORITHM 

A complete  listing  follows  of  the  simultaneous  iteration  algor- 
ithm for  generalized  symmetric  matrices  implemented  in  USA  Standard 
FORTRAN  as  SUBROUTINE  SIMITZ.  The  program  and  its  documentation 
are  separately  sequenced. 

■ 

I 
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SUBROUTINE  SIMITZ  (N,  P,  KM,  £PS,  IP,  OP,  INF,  EM,  X,  MN , 0,  WK> 


IDENTIFICATION 

SI  MIT  7 - ITERATIVE  COMPUTATION  OF  EIGENVALUES  LARGEST  IN  MAGNI- 
TUDE and  corresponding  eigenvectors  of  a real  general- 
ized symmetric  matrix 

fortran  SUBROUTINE  SUBPROGRAM 
US  AIR  FORCE  FLIGHT  DYNAMICS  LABORATORY 
MR  IGHT-PATT EPSON  AFB , OHIO  45433 
PURPOSE 

A REAL  N-SOUARE  MATRIX  C IS  O-S YMM E TP I C RELATIVE  TO  AN  N-SOUARE 
POSITIVE  DEFINITE  MATRIX  8 IN  CASE  9C  = C 'B  WHERE  C*  IS  THE 
TRANSPOSE  of  C.  GIVEN  AN  OPTIONAL  SET  OF  P INITIAL  APPROXIMATE 
EIGENVECTORS  OF  A REAL  N-SOUARE  O-SYMMETRIC  MATRIX  C CORRES- 
PONDING TO  P EIGENVALUES  OF  C LARGEST  IN  MAGNITUDE,  SIMITZ  COM- 
PUTES EM  EIGENVALUES  ANO  EM  CORRESPONDING  EIGENVECTORS  TO  A 
PRECISION  DEPENDENT  ON  THE  STRUCTURE  OF  8 ANO  C AND  ON  A GIVEN 
TOLERANCE  EPS.  THE  MATRIX  8 IS  PRESENTED  TO  SIMIT7  AS  AN  ALGO- 
RITHM FOR  CALCULATING  THE  STANOARO  INNER  PROOUCT  (W,  8?)  : W'BZ 
GIVEN  COLUMN  N-VECTORS  W and  Z IMPLEMENTED  AS  a FORTRAN  COM- 
PATIBLE REAL  FUNCTION  SUBPROGRAM.  THE  MATRIX  c IS  PRESENTED  AS 
A SUBROUTINE  SUBPROGRAM  WHICH  GIVEN  A COLUMN  N-VECTOR  Z CALCU- 
LATES ITS  IMAGE  w = CZ  UND-TR  THE  B-SYMMETRIC  MATPIX  C.  DEPEND- 
ING ON  THE  CHOICE  OF  B AND  C,  SIMITZ  APPLIES  TO  A WIDE  VARIETY 
OF  SYMMETRIC  eigenoroblems. 

CONTROL 

DIMENSION  X(MN,P>,  D(P>,  WKtK) 

INTEGER  P,  EM 
REAL  I® 

EXTERNAL  IP,  INF,  0° 


CALL  SIMITZ (N.  P,  X M , -PS,  IP,  OP,  INF,  EM,  X,  MN,  0,  W<» 

WHERE 

N IS  AN  INTEGER  INPUT  VARIABLE,  TH.  ORDER  Or  THE  MATRIX  C. 

P IS  AN  INTEGER  INPUT  VARIABLE,  THE  NUMBER  OF  SIMULTANEOUS 

ITERATION  VECTORS. 

KM  AS  AN  INTEGER  INPUT  VARIAPLE  IS  IN  MAGNITUDE  THE  MAXIMUM 

NUMBER  OF  ITERATION  STEPS  TO  BE  cXECUTEO.  IF  KM  IDENTIFIES 

A NEGATIVE  value  THEN  P INITIAL  APPROXIMATE  EIGENVECTORS 
ARE  ASSUMED  TO  BE  “RESENT  IN  THE  ARRAY  X.  OTHERWISE  SIMITZ 
SUPPLIES  RANDOM  INITIAL  EIGENVECTORS. 

KM  AS  AN  INTEGER  OUTPUT  VARIABLE  IDENTIFIES  THE  NUMBER  KS  OF 
ITERATION  STEPS  FINALLY  USEO  IN  THE  CALCULATION  BF  EM 
EIGENVECTORS. 

£R S IS  A REAL  INPUT  VARIABLE,  THE  TOLERANCE  FOR  ACCEPTING 
EIGENVECTORS. 

IP  IS  AN  EXTERNAL  INPUT  VARIABLE,  THE  N AM£  OF  A FORTRAN  COM- 
PATIBLE REAL  FUNCTION  SUBPPOGRAM  OF  TH£  FORM  IP(N,  Z,  W> 
WHICH  MUST  Frf(;pN  TH  - INNER  PROOUCT  (W,  37>  = W*17  OF  THE 
VECTORS  IDENTIFIED  BY  THE  N-ARPAVS  7 AND  W RELATIVE  TO  THE 
POSITIVE  DEFINITE  MATRIX  R. 

OP  IS  AN  EXTERNAL  INPUT  VARIABLE,  TH-  NAMl  Or  A FORTRAN  COM- 
PATIBLE SURPOUTINE  SUBPROGRAM  OF  THE  FORM  OP«N,  ',  W) 


SIMITZ 

z 

SIMIT7/0 

2 

SIMIT7/D 

3 

SIMITZ/C 

4 

SIMITZ/O 

5 

SIMITZZD 

6 

SIMITZ/D 

7 

SIMITZ/O 

8 

SIMITZ/D 

9 

SIMITZ/O 

10 

SIMITZ/O 

11 

SIMITZ/D 

12 

SIMITZ/O 

13 

SIMITZ/D 

14 

SIMITZ/O 

15 

SIMITZ/O 

16 

SIMITZ/D 

17 

SIMITZ/O 

18 

SIMITZ/O 

19 

SIMITZ/D 

20 

SIMITZ/O 

21 

SIMITZ/O 

22 

SIMITZ/D 

23 

SIMITZ/O 

24 

siMrrz/o 

25 

SIMITZ/O 

26 

SIMITZ/D 

27 

SIMITZ/O 

28 

SIMITZ/O 

29 

SIMITZ/D 

30 

SIMITZ/O 

31 

SIMITZ/O 

32 

SIMITZ/O 

33 

SIMITZ/O 

34 

SIMITZ/D 

35 

SIMITZ/D 

36 

SIMITZ/D 

37 

SIMITZ/O 

38 

SIMITZ/O 

35 

SIMITZ/O 

40 

SIMITZ/O 

41 

SIMITZ/O 

42 

SIMITZ/O 

43 

SIMITZ/O 

44 

SIMITZ/O 

45 

SIMITZ/D 

46 

SIMITZ/O 

47 

SIMITZ/O 

48 

SIMIT7/0 

49 

SIMITZ/O 

50 

SIMITZ/O 

51 

SIMITZ/O 

52 

SIMITZ/D 

53 

SIMITZ/O 

54 

SIMITZ/O 

55 

SIMITZ/O 

56 

SIMITZ/O 

57 
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WHICH  HIJST  CALCULATE  THE  IMAGE  W OF  TH£  VECT0®  1 1 £ NTIFI EO 
BY  THE  N-ARRAY  7 UNDER  THE  N-SQUAR  E MATRIY  C WITHOUT  OVER- 
WRITING 7. 

INF  IS  AN  EXTERNAL  INPUT  VARIABLE,  THE  NAME  OF  A FORr RAN  COM- 
°A  TIBLE  SURROJTINE  SUBPROGRAM  WHICH  MAY  BE  USEO  r OR 
OBTAINING  INFORMATION  OR  TO  EXERT  CONTROL  OUFING  ;X-.CUTION 
OF  SIMIT7.  INF  HAS  TH  F FORM  INF(KS,  G,  H,  F)  WHERE 

<S  IS  AN  INTEGER  OUTPUT  VARIABLE,  THE  NUMBER  Or  THE  NEXT 
ITEPATION  STE°. 

G IS  AN  INTEGER  OUTPUT  VARIABLE,  TH'  N'JMB.  R Or  ALR.ADY 
ACCEPTED  EIGENVECTORS. 

H IS  AN  INTEGER  OUTPUT  VARIABLE,  THE  NUMBER  Or  ALREADY 
ACCEPTED  EIGENVALUES. 

F IS  A REAL  OUTPUT  VARIABLE  P-ARRAY,  tRROF  QUANTITIES 
MEASURING  RESPECTIVELY  THE  STATE  OF  CONVERGE  NCE  OF 
THE  P SIMULTANEOUS  ITERATION  VECTORS,  IN  A10ITION, 

IF  CONVERGENCE  FAILS  IN  SUBROUTINE  IMTOL?  A r TiR  G 
EIGENVECTORS  HAVE  BEEN  ACCEPTED,  TH^N  FJG+ll  IS  RE- 
PLACED BY  UDO. *FLOAT(I£RR>  WHERE  IERR  IS  THE  ERROR 
INDICATOR  OUTPUT  BY  IMTOL2.  EACH  t LE  ME  NT  Or  THE  ARRAY 
F IS  INITIALLY  SET  BY  SIMTT7  TO  THE  VALUE  4.3. 

EM  AS  AN  INTEGER  INPUT  VARIABLE  IS  THE  NIIM8_ R OF  EIGENVALUES 
TO  BE  COMPUTED,  C ,LT.  EM  ■ LT,  P ,LE.  N . Li . MN. 

EM  AS  AN  INTEGER  OUTPUT  VARIABLE  IS  THE  NUMBER  OF  EIGENVECTORS 
COMPUTED  THROUGH  <M  ITERATION  STEPS, 

X AS  A REAL  N-BY-o  INPUT  ARRAY  IS  A Si T OF  P OPTIONAL  INITIAL 
APPROXIMATE  EIGENVECTORS  X(I,1>,  ....  X(I,P>,  1*1,  ..., 

N,  INTERPRETED  BY  SIMIT7  IF  KM  IS  NEGATIVE. 

X AS  A REAL  N-BY-d  OUTPUT  ARRAY  IS  A SET  OF  EM  EIS" N VtCTO RS 
X(I,1>,  ...,  X HI, EM) , 1 = 1,  ....  N,  COMPUTED  THROUGH 
IABS(KM)  ITERATION  STEPS  WITH  THE  REMAINDER  OF  X CONSISTING 
OF  P - £M  APPROXIMATE  EIGENVECTORS.  THE  N-BY-P  MATRIX  X 
SATISFIES  X'BX  = I,  THAT  IS,  THE  EIGENVECTORS  OR  C ARE 
3-ORTHOGONAL. 

MN  IS  AN  INTEGER  INPUT  VARIABLE  WHICH  IDENTIFIES  TH'  LEADING 
DIMENSION  IN  THE  CALLING  PROGRAM  OF  THE  ARRAY  X. 

0 IS  A REAL  OUTPUT  P-APRAY  OF  WHICH  0(1),  ....  O(.'M)  ARE  THE 
EIGENVALUES  OR  C LARGEST  IN  MAGNITUDE  IN  DE CR E AS TNG  ORDER 
CORRESPONDING  TO  THE  COMPUTED  EIGENVECTORS  X(I,1>,  ..., 

X ( I , EM)  , 1=1,  ....  N.  0 ( i M + l ) 0(0-1)  CONTAIN 

APPROXIMATIONS  TO  PROGRESSIVELY  SMALLE0  SUCH  E IG  rN V A LUES. 

D ( P)  CONTAINS  THE  MOST  RECENTLY  COMPUTED  VALUE  Oc  t,  WHERE 
THE  INTERVAL  (-L,  E)  IS  THE  INTERVAL  OVER  WHICH  ’Hf 
CHEBYSHEV  ACCELERATION  WAS  PERFORMED. 

WK  THE  INITIAL  LOCATION  OF  AT  LEAST  o»»2  ♦ 3»P  ♦ 2*N  = K 
CONSECUTIVE  STORAGE  LOCATIONS  WHICH  MAY  NOT  BE  OVER- 
WRITTEN WHILE  SIMIT7  IS  IN  EXECUTION. 

OTHER  PROGRAMMING  INFORMATION 

SI  MI T 7 EMPLOYS  a OATA  STATEMENT  TO  ASSIGN  TO  A MACHINE  DEPEN- 
DENT REAL  VARIABLE  MT  THE  QUOTIENT  OF  TH  _ SMALL-ST  POSITIVE 
REAL  VALUE  REPRESENTABLE  BY  FORTRAN  ANO  THE  SMALLEST  REAL  VALUE 
WHOSE  SUM  WITH  l.J  EXCEEDS  1.0. 

THE  PERFORMANCE  OF  SIMIT7  IS  STRONGLY  DEPEND. NT  UPON  THE  CHOICE 
OF  INPUT  PARAMETERS  ANO  UPON  THE  CAR.FUL  PREPARATION  OF  THE 
SUBPROGRAMS  IP  ANO  0®.  TME  USER  SHOULD  OEVELO®  EXP'OTENCE 
WITH  SIMIT7  ON  PPOBLEMS  OF  M00E  = ATE  SI7E  Oc.FORi  INVENTING 


SI  MI T 7/D  55 
SIMIT7/0  59 
SIMIT7/0  6 G 
SIMIT7/D  61 
SIMITI/O  62 
SIMIT7/0  63 
SIMIT7/0  6*, 
SIMIT7/D  65 
SIMIT7/0  66 
SIMIT7/D  67 
SIMITT/O  68 
SIMIT7/D  69 
SIMIT7/D  7j 
SIMIT2/0  71 
3IMIT7/D  72 
SIMIT7/0  73 
SIMITT/O  74 
SIMIT7/D  75 
SIMIT7/0  76 
SIMIT7/D  77 
SIMIT7/0  78 
SIMIT7/0  79 
SIMIT7/D  33 
SIMIT7/D  81 
SIMIT7/D  82 
SIMIT7/D  83 
SIMIT7/D  89 
SIMIT7/D  85 
SIWIT7/0  86 
SIMIT7/0  37 
SIMIT7/0  83 
SIMIT7/D  89 
SI  MI T 7 /D  93 
SIMIT7/D  91 
SIMIT7/D  92 
SIMIT7/D  93 
SIMIT7/D  99 
SIMIT7/D  95 
SIMIT7/D  96 
SIMIT7/0  97 
SIMIT7/D  98 
SIMIT7/D  99 
SIMI T7 /OlOL 
SIMIT7/D1Q1 
SIMIT7/D1I2 
3TMI T7/D1 3 3 
SIMIT7/0139 
SIMIT7/D105 
SIMIT7/01u6 
SIMIT7/0107 
SIMIT7/D1J8 
SIMIT7/D109 
SIMIT7/011. 
SIMIT7/0111 
SIMIT7/D112 
SIMIT7/DU3 
SIMIT7/0114 
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PROCESSOR  time  ON  VERY  LARGE  PROBLEMS  pop  WHICH  THE  ”=0Cl0UR2  SIMIT7/DU5 
IS  ULTIMATELY  INTENDED.  SIMIT7/D116 

OTHER  PROGPAMS  REQUIRED  STMIT7/0117 

FUNCTION  RANF  S I M I T 7/DI  l*i 

RETURNS  UNIFORMLY  0 I STJ IBUTr  0 PANOOM  NUMBERS  ON  THL  OPEN  SIMIT7/D119 
INTERVAL  (0,  1)  GIVEN  ANY  ONF  ARGUMENT  OF  ANY  TV":.  SIMIT7/012./ 

SUBROUTINE  TPE02  SIMIT7/D121 

IS  THE  EISPAC<  < <* » PROGRAM  WHICH  COMPUTES  A HOUSEHOLDER  SIMIT7/D122 

TRIDIAGONAL  FORM  OF  A REAL  SYMMiT°IC  MATRIX.  SIMIT7/D123 

SUBROUTINE  IMT0L2  SIMIT7/D12- 

IS  the  EIS°AC<  FROGRAM  WHICH  COMPUTES  the  EIGENVALUES  AND  EIMIT7/D125 
ORTHONORMAL  EIGENVECTORS  OF  A SYMMETRIC  TRIOIAGov'AL  MATRIX.  SIMIT7/D126 
FUNCTION  IP  SIMIT7/0127 

IS  DESCPIBEO  ABOVE.  SIMIT7/0A28 

SUBROUTINE  OP  SIMIT7/D12R 

IS  OESCRIREO  ABOVE.  SIMI  T 7 /m  3., 

SUBROUTINE  INF  SIMIJ7/0131 

IS  OESCPIBEO  ABOVE.  SIMIT7/D132 

METHOO  SIMIT7/DA  33 

SIMIT7  REPRESENTS  RESULTS  OF  EXTENSIVE  MODIFICATION*  AND  T-  STS  3IMIT7/D134 
OF  SUBROUTINE  RIT7IT  (1).  a'n  ANSI  F 0 s T RAN  TRANSLATION  OF  THE  SIMIT7/0135 

ALGOL  6C  PROCEDURE  Of  TH'  SAME  NAM  (31.  TH,  BASIC  p'ITISHAUSER  SIMIT7/0136 
-REINSCH  ALGORITHM  IS  RR: SERVED,  SIMIT7/0137 

REFERENCES  SIMIT77D133 

(1)  PAUL  J.  NIKOLAI  AND  NAI-KUAN  TSAO,  THE  A R L LIN  I'  ALGlBRA  SIMIT7/D135 

LIBRARY  HANDBOOK,  APL  TP  7*,-3106,  A;ROSRAC;  PIS' ARCH  (.ABO®-  STMIT7/T140 
ATOPIES,  WRIGHT-DATTERSON  AF  R,  OHIO,  1974.  3I“IT7/0141 

(2)  HEIN7  RUTISHAJSER,  COMPUTATIONAL  ASPECTS  0r  F.  L.  BAU.R'S  SIMIT7/0142 

SIMULTANEOUS  ITERATION  METHOD,  NU"  ER . MATH.  13(1  *><>9),  4-13,  SIMIT7/n:43 

(3)  , STMULTAN-.OUS  ITERATION  M-THOO  pos  SYM-  SIMIT7/D144 

METRIC  MATRICES,  NUMtR.  MATH,  16(19701  , 236-223.  SIMIT7/0145 


(41  B.T 

. SMITH  et  al, 

MATRIX  EIGLNSYST  M 

POUT  I U 

_ S-  .1 "OftCK 

SI  MI T 7 /P 

Nh 

GUIDE,  Lr.CTU°E  NOT 

ES  IN  COMPUT.R  SCI 

NCE  5, 

S°-  I NGiR- 

V .'LAG 

SIMIT7/D 

147 

NEW 

YORK,  1974. 

5IMIT7 /0 146 
SIMIT7/P149 

EXTERNAL 

INF,  IP, 

OP 

SIMIT7 

4 

INTEGER 

EM,  G, 

H,  I . 

IG, 

IK  . 

J 9 

S I M I T 7 

c 

* 

JK,  JO, 

K,  KM, 

KS, 

L , 

LF, 

SIMIT7 

* 

* 

LI,  M, 

MN,  Ml, 

N, 

o f 

71, 

SIMIT  7 

7 

* 

72 

SIMITT 

? 

LOGICAL 

ORIG 

5 I M I T 7 

9 

REAL 

D,  E, 

EPS,  El, 

_ 2 . 

ip  . 

MT, 

SIMIT7 

1C 

* 

S,  T, 

WK,  X 

SI  MIT  7 

11 

DIMENSION 

X (MN ,11  , n(  1)  . 

WK  (O  ,0 , 1> 

SI  MI T 7 

12 

S I M I T 7 

13 

OATA  MT  /. 

2203606416450  62 

£.-?7P/ 

SIMIT 7 

14 

THE  LOCAL  VAPIABLE  ARRAY*  FROM  (3)  AO  ASSIGN. 0 TO  IN- 
VARIABLE A R o AY  WK  AS  FOLLOWSI 


WXII.J.II  = 0(1, j)  , I,  J = l 
WK  (I,  1, 2)  = CX(  I)  , I = i.  , . 
WK( I , 2, 2)  = F ( I ) , I = 1,  ... 
WK( I, 3,21  = RQ ( I I , I = 1 , . . 
WK  (1,4,2)  s II  (II,  I = 1,  ... 
WK(I+N,4,2)  : Hill,  I ! 1,  , 
NOT  NEEOEO  ARE  V(I)  , I = 1 . 


R(  I)  , I o 


.,,,  P,  A NO 


S I M I T 7 

SI  M IT  7 
SIMIT7 
S I M I T 7 

SI  MI T 7 
SIMIT7 
SIMIT7 
S I M I T 7 
SI  MI T 7 
SIMIT7 
STM  I T 7 
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01 1 * J)  • I • J - 1 * •••• 

SIM  IT  7 

26 

SI  MI T 7 

27 

THE  NEXT  STATEMENT  15  STAPT. 

SI  MI T 7 

28 

£ = .0?*30 

SIMIT7 

29 

0=0 

SIMITZ 

30 

10  = 1 

SI  MI T 7 

31 

H = 0 

SI  MI T 7 

32 

= 1 i a 

SIMITZ 

33 

72  r 0 

SIMITZ 

34 

KS  = 0 

SI  MI T 7 

35 

1=1 

SIMITZ 

3b 

MK(P,1,2>  = .0E+00 

SIMITZ 

37 

00  10  L = 1,  0 

SIMITZ 

38 

WK(L,2,2)  = .4E+C1 

SIMIT7 

39 

WK(L,3,2>  = .C-+G0 

SIMITZ 

4 C 

10 

CONTI NUE 

SIMITZ 

41 

IP  (<M>  50,  50,  20 

SIMITZ 

42 

20 

90  43  L * It  » 

SIMITZ 

43 

DO  3C  1 = 1,  N 

SI  MI T 7 

44 

X(J,L)  = « 2‘  ♦Cl*9ANr(.  2_*C1>  - 

,1.*l1 

SIMITZ 

45 

30 

CONTINUE 

SIMITZ 

46 

4 0 

CONTINUE 

SIMITZ 

47 

50 

KM  r IA9S(KM) 

SI  MI T 7 

48 

ASSIGN  63  TO  IK 

SIMITZ 

49 

UP  = IG 

SIMITZ 

6 j 

LI  r o 

SIMITZ 

51 

GO  TO  990 

SIMITZ 

52 

PA TLE ICH-9I T7  STE® 

SIMITZ 

53 

STATEMENT  60  IS  LOOP. 

SIMITZ 

54 

6G 

DO  80  K = IG,  P 

SIMITZ 

55 

CALL  0°(N,  X(1,K>,  WK(1, 4.2)1 

SIMITZ 

56 

DO  70  J » 1,  N 

SIMIT7 

57 

X(J,K>  = WK( J,4,2) 

SIMITZ 

58 

70 

CONT  INIJE 

SI  MIT? 

59 

• 0 

CONTI  NUE 

SIMITZ 

6* 

ASSIGN  93  TO  IK 

SIMITZ 

61 

LP  = IG 

SIMITZ 

62 

LI  = ° 

SIMITZ 

63 

GO  TO  990 

SIMIT7 

64 

90 

IP  (KS)  150,  1GI,  150 

SIMITZ 

65 

MEASURES  AGAINST  JNHAPPY  CHOIC- 

OF  INITIAL  VlCTOPS 

SIMITZ 

6b 

100 

00  133  K = 1,  P 

SIMITZ 

67 

IF  <WK<K,K,i»)  173,  110,  132 

SIMIT7 

68 

110 

DO  120  1=1,  N 

SIMITZ 

69 

7(1, K)  = ,?:+01*9ANP(,2:  + 01)  - 

,1i*C1 

SIMITZ 

7* 

123 

CONTINUE 

SIMIT7 

71 

KS  = 1 

SI  MI T 7 

72 

130 

CONTI NUf 

SIMITZ 

73 

IP  (KS  - 1)  150,  140,  1^0 

SIMITZ 

74 

140 

ASSIGN  60  TO  IK 

SIMITZ 

75 

LP  = 1 

SIMITZ 

7b 

LI  = 0 

SIMITZ 

77 

GO  TO  993 

SIMITZ 

78 

150 

00  180  K = IG,  P 

SIMITZ 

79 

DO  170  L = K,  P 

SIMITZ 

9b 

S = ,0£*00 

SIMITZ 

81 

00  160  I = L,  P 

SIMITZ 

82 

23 
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S = S ♦ 44( 1,4,  1) *44  (I  ,L  ,1) 

SI H I T 7 

83 

160 

CONTINUE 

SIMIT7 

94 

44(L.K.l)  = -9 

SIHIT7 

95 

17j 

CONTINUE 

SIMIT7 

96 

H3 

conti  wr 

SIHIT7 

97 

call  T9ED2(P,  p - f.,  44(19, IS, 11,  0(191,  44(1,4,2),  44  ( 

TG.IG.l)  1 

SI  HI T 7 

98 

CALL  IMT9L2(o,  p - G,  0(191,  44(1,4,21,  4 4 ( If,,  IG , 11  , U 

SIHIT7 

89 

44(19,2,21  = AMAX1(44(IG,2,2)  , .l£*04*PLOaT  (LI  1 

3IHIT7 

90 

oo  no  < = ig,  p 

SIHITZ 

91 

n (< t = sort (aua«K-n(Ki,  .or *00 » i 

SIHIT7 

92 

190 

CONTINUE 

SIHIT7 

93 

REORDERING  EIGENVALUES  aNO  £ IG : NVEC  TOR  S ACCORDING  TO 

S I 7 » OF 

SIHIT7 

94 

THE  FORMER  II  ACCD“PLISM.‘0  IN  SUBROUTINE  IHTOL2. 

CIHIT7 

95 

DO  210  J = 1,  N 

SIHITZ 

96 

OO  210  4 = IG.  P 

SIHITZ 

97 

S = ,JE*3) 

SIHIT7 

98 

OO  2)0  L = IG,  p 

SIHITZ 

99 

S = S ♦ *(J,ll»44(L,4,l) 

SIHITZ 

1 j. 

?00 

CONTINUE 

SIHIT7 

101 

44(4,4,2)  = 9 

SI  HI T 7 

102 

’10 

CONT IN JE 

SIHITZ 

133 

OC  22^  4 = IG,  D 

SIHIT7 

104 

* (J.4)  = 44(4,4, 21 

SIHITZ 

1C5 

2 20 

CONT INUE 

SIHIT7 

106 

230 

CONTINUE 

si-itz 

1JZ 

4S  = <S  ♦ 1 

SIHITZ 

1 3 8 

£ r AM  A*1 (0(“ 1 , E 1 

SIHITZ 

109 

®AN00HI7ATT0N 

SIHITZ 

110 

IF  (3  - 7ii  260,  240,  240 

SIHITZ 

111 

240 

OO  25  0 J = 1,  N 

SIHITZ 

112 

X(J,P1  = ,2- ♦C1*94NF(.2£*C1)  - .17+01 

SIHITZ 

113 

29) 

CONTINUE 

SIHITZ 

1 1 4 

JP  r D . l 

SIHITZ 

115 

ASSIGN  260  TO  14 

SIHITZ 

116 

LF  s P 

SIHITZ 

117 

LI  i ° 

SIHITZ 

113 

GO  TO  990 

SI  HI T 7 

119 

COMPUTE  CONTROL  QUANTITIES  04(11. 

SIHITZ 

12C 

260 

OO  313  4 = IG.  JP 

SIHITZ 

121 

S = (0(41  - t)*(P(4l  ♦ El 

SIHITZ 

122 

IF  (Si  270,  270,  2*0 

S IH I T 7 

123 

27C 

44(4,1,2)  * « » E *0C 

SI  HI T 7 

124 

GO  TO  310 

S I H I T 7 

125 

2 90 

IF  ( E 1 3 CO , 290  , 30  0 

SIHITZ 

126 

293 

44(4,1,2)  = • 1 E ♦ 0 4 ♦ ALOG  (0(4)) 

SIHITZ 

127 

GO  TO  310 

SIHITZ 

128 

300 

44(4,1,2)  - A L 0 G ( ( O ( 4 1 ♦ SOPT (S 1 ) /E ) 

SIHITZ 

129 

31) 

CONTI NU£ 

SIHITZ 

13b 

ACCEPTANCE  TEST  FOR  IGE'IVALUES  INCLUDING  ADJUSTMENT 

OF  iM  ANO 

SIHITZ 

131 

M SUCH  THAT  O ( EM|  ,GT.  t,  O(H)  ,GT,  _ AND  Q(_M)  OOi  S 

NOT 

SIHITZ 

132 

OSCILLATE  STFONGLT 

SIHITZ 

133 

I * 71  - 1 

SIHITZ 

134 

K * tt 

SIHITZ 

135 

*23 

4 * 4 ♦ 1 

SIHITZ 

136 

IF  (£M  - 4)  373,  330,  330 

SIHITZ 

137 

330 

IF  (0(41  - -.)  36C,  362,  34C 

SIHITZ 

133 

340 

IF  (I)  320,  320,  350 

SIHITZ 

139 

II 
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353 

IF  (0(0  - .999E*0;*WK(<,3,2>)  *6C  , 360  , 32C 

SIMIT7 

1.: 

360 

CONTI  NUF 

SIMIT7 

141 

EM 

1 * K - 1 

SI  MI T 7 

142 

STATEMENT  370  tS  -:X4. 

SIMIT7 

143 

370 

15 

(EM)  380  , 1130  , 380 

SI  MI T ? 

194 

380 

K 

= H 

SI  MI T 7 

145 

S 

= , IE ♦ 0 1 ♦ ,lfc+C0*ERS 

SIMIT7 

146 

390 

K 

= K * 1 

SIMIT7 

147 

IF  (OIK))  4 C 0 , .10,  4C3 

SI  MI T 7 

148 

<•00 

IF  (0(0  - S»HK(K,3,2>>  39Q,  390,  410 

SI  MI T 7 

149 

<♦10 

CONTI NUE 

SIMIT7 

150 

H 

= K - 1 

3IMIT7 

151 

K 

= EM 

SI  MI  T 7 

152 

«*20 

< 

= K ♦ 1 

SIMIT7 

153 

IF  (K  - H)  430,  433  , 450 

SIMIT7 

154 

430 

IF  (0(0  * E)  (,(,£,  440  , 4 2C 

SIMIT7 

155 

443 

CONTI NUE 

SI  MI T 7 

156 

H 

= < - 1 

SIMIT7 

157 

ACCEPTANCE  TEST  FOR  EIGENVECTORS 

SIMIT7 

153 

450 

L 

= G 

SIMIT7 

159 

E2 

! = • 0 E ♦ 0 0 

SIMIT7 

16: 

00  590  < = IG,  jo 

S I M I T 7 

161 

IF  ( K - ( L ♦ 1))  51 C , 460,  510 

SIM  I T 7 

162 

CHECK  FOR  NcSTEO  EIGENVALUES 

SIMIT7 

163 

463 

L 1 K 

SIMIT7 

16.. 

11  = K 

SIMIT7 

165 

S = .5E+00/FLOAT (KS) 

SIMIT7 

166 

T = .1E*01/FLOAT(KS»M) 

SIMIT7 

167 

470 

L = L ♦ 1 

SIMIT7 

163 

IF  (L  - J°)  480,  480,  490 

S I M I T 7 

169 

480 

IF  (HK(L, 1 ,?> * (H< (L , 1, 2)  ♦ S)  ♦ T - HO  L- 1 , 1 , 2 > *HK  ( L- 1 . 1 ,2  ) > 

SIMIT7 

173 

* 

490,  49C , 470 

SIMIT7 

171 

490 

CONTINUE 

SIMIT* 

172 

L = L - 1 

SIMIT7 

173 

THE  NEXT  STATEMENT  IS  E X 5, 

S I M I T 7 

174 

IF  (L  - HI  5 1C , 510,  500 

EIMIT7 

175 

503 

L = LI  - 1 

SIMIT7 

176 

GO  TO  600 

SIMIT7 

177 

510 

CALL  0 3 ( N , X(1,KI,  HOI, 4, 2)) 

SI  MI T 7 

178 

S = • u£*0C 

SIMIT7 

179 

00  54 C J * 1,  L 

STMIT7 

181 

IF  ( A9S  (0  ( J>  - 0(0  ) - .1,-01*0101  520,  54,,  54' 

SI  MIT  7 

181 

520 

T = I R ( N,  WK (1,4,2) , X ( 1 , J ) ) 

SIMIT7 

182 

00  530  I = 1,  N 

SIMIT7 

183 

HOI, 4, 2)  = WO  I , . ,2  > - T*X  ( I , J) 

SIMIT7 

13. 

530 

CON'*  INUE 

SIMIT7 

185 

S = S ♦ T*T 

SIMIT7 

136 

540 

CONTINUE 

SIMIT7 

187 

T - . E * 3 1 

SI  MI T 7 

1 S3 

IF  (S  ,N£.  . OE  * 00 ) T = I R ( N , W*<1,4,2I,  H<(1,4,2)1 

SIMIT7 

189 

E2  = AMAXKE2,  SORT(T/(S  ♦ T * ) > 

SIMIT7 

193 

IF  (K  - L)  590,  550  , 590 

SIMIT7 

191 

TEST  50R  ACCEPTANCE  OR  GROUP  OF  EIGENVECTORS 

SIMIT7 

192 

550 

IF  (L  ,G£.  EM  ,AN9,  0(EM)  *WO£M,2,2>  ,LT.  ERS»(0(£M>  - Ell 

SIMIT7 

191 

• 

G = £M 

SIMIT7 

194 

IF  (£2  - MK ( L , 2 , 2) ) 563,  58,,  583 

SIM  I T 7 

195 

560 

00  573  J = LI,  L 

SIMIT7 

196 

25 
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WK(J.2,2>  = 3 2 

SI9IT7 

197 

^7; 

CONTINUE 

SI9IT7 

198 

5 83 

IF  (L  .LE.  3 M , AND.  OIL)  (L,  2.2)  ,LT.  cPS»(0(L)  - )>  G = L 

SI9IT7 

199 

593 

CONTINUE. 

SI9ITZ 

200 

ADJUST  M, 

SI9IT7 

201 

STATEMENT  633  IS  *6. 

SI9IT7 

202 

8 0 3 

IS  = S ♦ 1 

SI  91 T 7 

203 

Ic  ( £ - • 4i«0 1*9(11)  61„,  610,  623 

SI9IT7 

204 

813 

9=1 

SI9ITZ 

205 

K = 1 

SI9ITZ 

236 

GO  TO  633 

SI9IT7 

207 

623 

E2  = .2'iOl/E 

SI9IT7 

203 

51  = .51I+00*l2 

SI  MIT  7 

239 

< = 2*  IN  T (,4£*31/AMIN1 (MKCl.1,2)  . .4t>Cl>) 

SI  91 T 7 

210 

9 = MIN.  (9,  K) 

SI  91 T 7 

211 

REDUCE  £9  IF  CONVERGENCE  WOULD  9£  TOO  SLOW. 

SIMIT7 

212 

6 33 

Ic  (W <(.9,2,21)  64C,  693.  640 

SI9IT7 

213 

643 

IF  (FLOAT  ( Rf  ) - .9v*OC»FLOAT(KM)  ) 653  , 690  , 69|. 

SI  9 1 T 7 

214 

650 

S = FLOAT(K)»WK(E-,l,2) 

SI9IT7 

215 

IF  (S  - .5E-01)  663 , 67C , 673 

SI9IT7 

216 

663 

T = ,5E+30*S*WK(F9,1,2) 

SI9IT7 

217 

GO  TO  683 

SI  91 T 7 

218 

673 

T = W<(£9,1,2>  ♦ ALOS  (•  5Z*  0 0 ♦ . 5E*U0  * Z.  X“  ( -.  2c +C  1*S>  > /F  LOAT  (K) 

SI9IT7 

219 

683 

S = AL0G(U(Z9>»WK(FM,2,2I/(;PS»(0(£9>  - £))) 

SI9IT7 

220 

IF  (S»fLOAT(i<S)  ,r,T.  T*FLOAT((K9  - KS)*K9)>  E9  : ;9  . 1 

SI9IT7 

221 

STATEMENT  690  IS  :X2. 

SI9IT7 

222 

693 

OO  733  K = IG,  JP 

SI9IT7 

223 

WK(<,  9,2)  = 0(0 

SI9IT7 

224 

703 

CONTI NU£ 

SI9IT7 

225 

CALL  INF ( KS , G,  H,  W<(1,»,JII 

SI9IT7 

226 

IF  (G  .or.  t9  .09.  <5  . G£.  KM)  GO  TO  1133 

SI9IT7 

227 

STAT-.9-ENT  713  IS  f«l. 

SI9IT7 

228 

713 

IF  (<S  ♦ 9 - KM)  73C , 733  , 723 

SI  91 T 7 

229 

720 

7.2  = -1 

SI  91 T 7 

230 

IF  (9  , GT , 1)  9 = 2»(<KU  - KS  ♦ 1)7?) 

SI9IT7 

231 

7?a 

Ml  = 9 

SI  91 T 7 

232 

S90PTCUT  LAST  INTO  » M _ OI  A T£  9L0CK  IF  ALL  F(I)  A c E SUFFICIENTLY 

SI9IT7 

233 

SMALL. 

SI9IT7 

2 34 

Ic  (L  - £9)  73  L , 7 43 , 74  3 

SI9IT7 

235 

743 

S = 0«F9)*WK(£M,2,2l/(  PSMI3I  .9)  - •_)) 

SI  91 T 7 

236 

T = S*S  - .IE* Cl 

SI9IT7 

237 

IF  IT)  63  . 60  , 760 

SI9IT7 

236 

750 

S = ALOGIS  ♦ SOP  T (T)) / ( WK ( £ 9 , 1,2)  - W<(H+1, 1,2)1 

SI  91 T 7 

239 

91  = ?*INT(.5:.*C0»S  ♦ .1J1£‘01) 

SI9IT7 

24. 

IF  (91  - 9)  770,  77r , 763 

SI  91 T 7 

241 

760 

91  = 9 

SI  91 T 7 

242 

GO  TO  763 

SI9IT7 

243 

770 

72  = -1 

SI9IT7 

244 

CM  :OY>9LV  TTc’ATION 

SI  91 T 7 

245 

780 

IF  (m  - 1)  S0G,  700.  620 

S 1 9 1 T 7 

246 

793 

00  81 C K = IG,  ° 

SI  91 T 7 

247 

CALL  O’lN,  X(j,K),  WK 1 1, 4 ,2) ) 

SI9IT7 

248 

00803  I s 1 , N 

SI  91 T 7 

2 49 

X(I,K)  = WK(I,4,2) 

SI9IT7 

25  j 

803 

CONTINUE 

SI9IT7 

251 

*IC 

CONTINU-. 

SI9IT7 

252 

GO  TO  9l  3 

SIM  I T 7 

253 

o o o 
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820  LI  = Ml  - 4 

DO  693  K * IG,  p 

CALL  0°<N,  XII, K),  WK<1,4,2I) 

00  83C  I = li  N 

IK  = I + N 

WK(K,4,2»  = E1**K(I,4,2J 
830  CONTINUE 

CALL  0°(N,  HK(N*1,4,2>  , WK(», <,,?)) 

00  840  I = 1,  N 

X(I,K>  = E 2* WK ( 1 , 4 « 2 ) - X(I,K) 

840  CONTINUE 

IF  <L1>  89C,  "50,  85  0 
850  00  880  J = 4,  Ml,  l 

CALL  0»(N,  X(1,K),  WK(1,4,2>> 

DO  "60  I = 1,  N 

IK  = I * M 

MKIIK.4,21  = EZ»WK<I,4,2)  - WKlIK.u,?! 
860  CONTINUE 

CALL  OP(N,  WK(N+1,4,2I,  WK(1,4,2>> 

00  870  I * 1,  N 

X(I,K)  r S2*WK<I,4,2>  - X(I,K) 

870  CONTINUE 

880  CONTINUE 

890  CONTINUE 
900  ACS IGN  910  TO  IK 
LF  = IG 
LI  - ° 

GO  TO  993 

C DISCOUNTING  THE  ERRO  = OUANTITIlS  F 


910 

IF  (G  - H) 

92.,  97C  , 

97 1 

923 

IF  (M  - 1) 

950  , 930, 

950 

930 

00  943  K = 

IG,  h 

WK(  K , 2 , 2) 

R WK(K,2 

,21*  ID  (H*  n/O  (KM 

940 

CONTINUE 
GO  TO  973 

950 

T = EXP(-FLOAT t«l ) »W< ( H*l, 1,?) ) 

00  96  0 < = IG,  H 

s = ex°c-floatcmi>*(v4k(k,i.7>  - MK(H«i,i,?»n 

WK(K,2,2I  = S*NK«K,2,2>»l.l:*01  ♦ T»T>/l.lF»Ci  ♦ IS*T**<S*T>) 
960  CONTINUE 
970  KS  = KS  ♦ Ml 
Z2  = 72  - Ml 

C P0SSI9LE  REPETITION  ne  I MTt RM£OI AT E "TtPS 

IF  (7  2 > 980,  T1C,  713 
980  71  = 71  ♦ 1 

72  = 2*71 
M = M ♦ M 
GO  TO  6C 

PERFORMS  OP.TMONORMALI7ATION  OF  COLUMNS  1 THROUGH  LI  TF  ARRAY 
X ASSUMING  THAT  COLUMNS  1 THPOUOH  Lc  - 1 APE  ALR-AOY  ORTHO- 
NORMAL 


993 

00  1120  K = Lc , 

LI 

ORIG  = .TRUE, 

1 GO  0 

T = .:E*3t 

JK  R K - 1 

IF  ( JK 1 1043, 

1 C 40  , 

1010 

DO  1033  1=1 

. JK 

s : m : t 7 

25-* 

SI M 1 1 7 

7 55 

S I M I T 7 

25b 

SI  MI T 7 

257 

SIMJT7 

258 

SI  M I T 7 

259 

S I M I T 7 

26. 

S I M I T 7 

2 51 

SI  MI T 7 

262 

SIMIT7 

26? 

SIMIT7 

264 

SIMIT7 

2 65 

SIMIT7 

766 

SI  MI T 7 

267 

SI M I T 7 

268 

S I M I T 7 

269 

SI  MIT  7 

27>. 

SI  MI T 7 

771 

SIMIT7 

772 

SIMIT7 

273 

SIMIT7 

274 

S I M I T 7 

27s 

SI  Ml T 7 

776 

SI MI T 7 

777 

S I M I T 7 

7 7 

SIMIT7 

779 

SI  MI T 7 

78. 

SIMIT7 

281 

S I M I T 7 

782 

SIMIT7 

285 

SI'*IT7 

784 

SIMIT7 

7 85 

S I M I T 7 

786 

SI  MI T 7 

2 6 f 

S I M I T 7 

7 88 

S IMI T 7 

7 8 5 

SIMIT7 

£.9'. 

SI  MI T 7 

291 

SI  HI  T 7 

292 

SI  MI  T 7 

-93 

ST  M I T 7 

7 94 

SI  MI  T 7 

2 95 

S I h I T 7 

296 

SIMIT7 

297 

SIMIT7 

7 98 

SI M I T 7 

2 99 

SIMIT7 

33- 

SI MI T 7 

3J1 

SIMIT7 

3 2 

SI “I T 7 

iu  3 

SUIT* 

334 

SI  MI T 7 

305 

SI  HIT  7 

3 j6 

SI  MI T 7 

327 

SIMITZ 

338 

S I M I 7 7 

33  9 

SIMIT7 

31 . 

27 
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102C 

101C 

104G 


1054 
1 0 SO 

1070 

1090 


1090 

not 

1110 

1130 

C 

1 1 1C 

c 


1 14  C 


1150 

1150 


117C 

HOC 


1190 
120  0 
1210 


s = ip(n,  v « i ,i),  7(1, <n 
ip  (o=ir.)  mk(x,i,ii  = 5 
r = t ♦ s»s 

99  1 0 20  J = 1,  N 

X(J,X>  = v<  J,X>  - S»X(J,I) 
CONTINUE 
CONT  IN'Il 

S = IP(N,  X ( 1,7)  , X I 1 ,<)  ) 

T = S ♦ T 


IF 

(S  - ,1E-01*T) 

io*>;  • 

1.5.,  1050 

IF 

(T  - MT)  i:*3. 

1 c * 

i „ * ; 

ORI 

G = .FALSE. 

IF 

(S  - MT)  l?7u. 

1 u 7l  t 

i j: 

S = • b E ♦ 30 
S = SO*  T ( 5 1 
WX(X,x,  l>  = $ 

15  (SI  1090.  llOi,  1090 
S - . 1E+01/5 
DO  111;  J = 1,  N 
X( J,X)  = S * x ( J,K) 

CONTINUE 
CONTI  NIC 

CO  TO  IX,  (50.  90,  250,  910,  114;> 

STATEMENT  11X0  IS  EX, 

i"  - 5 

SOLV--  EICENVf-LU"  ®’OPL.-M  Oc  P=O.UCTION  Q5  MATRIX  C. 

ASSIGN  1 1 4C  TO  IK 
L5  = 1 
LI  = J° 

GO  TO  990 
OO  1154  X = 1,  .IP 

CALL  0 ° ( N , X (i ,<) , X (1 ,*)) 

00  1150  I = 1,  X 

«<(X,I,1)  = * I P ( N , X (1  , T(  , X(1,P|) 

CONTINUE 
CONTI  NUE 

CALL  TRi02(P,  JP,  MX,  0,  M<(l,4,?l,  MX) 

CALL  IMTOL?(P,  J*.  0,  MX (1,4, 21,  MX,  L) 

MX ( IG , 2 . 2 ) = ANfXl(MX(TG,?,2l  , . 1 E* 14* FLOAT ( L I ) 

ARRANGE  EIGENVALUES  IN  OPOtP  0*  OECR^ASING  ABSOLUT:;  VALUc. 
00  1213  1=1,  J® 

X = J 

00  1171  I = J,  JP 

15  ( AOS (0 ( I > ) ,GT.  A 9S ( P ( X ) ) ) K = I 
CONTINUE 

IF  (X  - J)  12C0 , 12CC.  11»C 
T = O(X) 

n(<)  = oi  ji 
0(J)  = T 

00  1190  1 = 1,  JP 
T = MX  ( I,  X , 1 ) 

MX(I,X,1)  = MXII.J.l) 

MK(I,J,1>  = T 
CONTINUE 
0(J>  = -0(.J) 

CONTI  NU: 

00  125C  J = 1,  N 


SIMIT7 
SI  HI T 2 
SIMIT7 
SIMITZ 
SIMITZ 
SIMIT  7 
SIM1T7 
SIMITZ 
SI-ITZ 
SIMITZ 
SIMITZ 
SI-ITZ 
SIM  I T 7 
SIMITZ 
SIMIT7 
SIMITZ 
SIMITZ 
3IMIT7 
SIMIT7 
SIMIT7 
SIMITZ 
SIMIT7 
S I M I T 7 
SI  m I T 7 
SI  MI T 7 
SIMIT7 
SI  MI T 7 
SIMIT7 
SIMIT7 
SIMIT7 
St  MI T 7 
SI  MI T 7 
SIMIT7 
SIMIT7 
SIMIT7 
SIMIT7 
SIMITZ 
SIMITZ 
SIMITZ 
SIMIT7 
SIMITZ 
SIMITZ 
SIMITZ 
SIM  I T 7 
SIMITZ 
SIMITZ 
SIMITZ 
SIMIT7 
SIMITZ 
SIMIT7 
SIMITZ 
SIMITZ 
SIMITZ 
SIMITZ 
SIMIT7 
SIMITZ 
SIMITZ 
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00  12  7 0 1=1.  1 0 

s = 

00  1220  K = 1 i J3 

S = S ♦ XI J,K)»MK(K, 1.1) 
122C  CONTINUE 

UK  < I «<*  1 2)  = S 
1230  CONTINUE 

DO  121.3  I - 1,  JP 
X (J,I>  = NK(I,U,21 
1240  CONTINUE 
1250  CONTINUE 
KM  = KS 
O(P)  = E 
RETURN 
ENn 


SIMITZ 

365 

SI  MI T 7 

369 

SIMITZ 

37  u 

S I M I T 7 

771 

ST  M I T 7 

372 

S I M I T 7 

373 

SIM  IT  7 

374 

SIMITZ 

775 

SI  MI T 7 

3 76 

SI  MIT  7 

377 

SIMITZ 

376 

SIMITZ 

379 

SI mI T 7 

33  j 

3 1 M 1 T 7 

381 

SI  MI T 7 

732 

29 
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